 /**************************************************************************//**
 * @file    Siklon.c
 * @author  İsmail SEZEN (email: sezenismail@gmail.com)
 * @date    06.02.2011
 *
 * <b>Description</b>
 *
 * Most of the work done here. This module contains routines to calculate
 * the cyclone centers located in the charts represented by netCDF data format.\n
 *
 * Copyright (C) 2009 İsmail SEZEN
 *
 * <b>License</b>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 *****************************************************************************/

/* Include Directives */ //{
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <netcdf.h>
#include "Helpers.h"
#include "CycloneDoubleLinkedList.h"
#include <array_helper.h>
#include <logger.h>
#include "AttributeHelper.h"
#include "siklon_internal.h"
#include "NC_Helper.h"
//}

/* Define Directives */ //{

/*! \brief Elimination radius in Km.
 *
 * this constant is used to eliminate the shallow cyclones around the deeper one
 */
#define ELIMINATION_RADIUS 300

//}

/* Macros */ //{

/*! \brief looping helper macro \a i is loop var \a x is start value and \a y is decision variable */
#define loop(i,x,y) for( (i)=(x); (i)<=(y); (i)++ )
/*! \brief returns the length of an array */
#define length(a) (sizeof a / sizeof a[0])
/*! \brief swaps \a a and \a b values*/
#define SWAP(a, b) {double t=a;a=b;b=t;}

//}

/* Private Routines ---------------------------------------------------------------------------------*///{

/*! \var cyc
 *  \brief Holds all cyclones that found
 *
 * This variable holds the cyclones that found in the level for a time step.
 * For every time step this variable is reset and reuse.
 */
struct CYCLONE cyc;

 /*!
 * \brief Connects two double linked list
 *
 * this function is used to connect two double linked list contains found cyclones.
 *
 * \param[in] dl1_first first element of double linked list 1
 * \param[in] dl1_last last element of double linked list 1
 * \param[in] dl2_first first element of double linked list 2
 * \param[in] dl2_last last element of double linked list 2
 * \return void
 * \note this is a private helper function
 * \warning note that parameters are pointer of pointer
 */
static void connect_dls(struct cyclone_item** dl1_first, struct cyclone_item** dl1_last,
						struct cyclone_item** dl2_first, struct cyclone_item** dl2_last ) {

	if(*dl1_last) {
		(*dl1_last)->next = (*dl2_first);
		(*dl2_first)->prev = (*dl1_last);
		(*dl1_last) = (*dl2_last);
	}else{
		(*dl1_first) = (*dl2_first);
		(*dl1_last) = (*dl2_last);
	}
}

/*!
 * \brief Cheks the given initial time parameters against NC file's time range
 *
 * this function is used to check validity of user-given initial time index parameters and
 * compare them against \a info.dims.time.len value.
 *
 * \param[in] start this is start time given by the user.
 * \param[in] end This is end time index given by the user.
 * \param[in] time_step his is time step given by the user.
 * \return void
 * \note this is a private helper function
 * \warning note that parameters are pointer
 */
static void checkTime(int* start, int* end, int* time_step) {
	int len = info.dims.time.len-1;
	if(*start<0 || *start > len ) *start = 0;
	if(*end<0 || *end > len ) *end = len;
	if(*start > *end) SWAP(*start, *end); /*Swap values*/
	if(*time_step < 1 || *time_step > len ) *time_step = 1;
}

/*!
 * \brief Cheks the given initial latitude/longtitude start and end values against files's limits
 *
 * this function is used to check validity of initial latitude/longtitude parameters and
 * compare them against NC file's domain limits. if given values are ou of range, they
 * have been automatically set to \a range.min and \a range.max values.
 *
 * \param[in] start This is start value of latitude/longtitude given by the user.
 * \param[in] end This is end value of latitude/longtitude given by the user.
 * \param[in] range This is valid latitude/longtitude range handled from the NC file.
 * \return void
 * \note this is a private helper function
 * \warning note that parameters are pointer
 */
static void checkLatLon(float* start, float* end, struct range_float range) {
	float min = range.min;
	float max = range.max;
	if(*start<min || *start > max) *start = min;
	if(*end<min || *end > max) *end = max;
	if(*start == *end) {
		*start = min;
		*end = max;
	}
	if(*start > *end) {//Swap values
		SWAP(*start, *end);
	}
}

/*!
 * \brief Sets the values of internal domain variables by initial variables entered by the user
 *
 * This function is used to set the internal variables to be used by a variaety of functions in the project.
 * Mainly, initial latitude/longtitude/time variables is set here according to user input.
 *
 * \param[in] lat_start This is start value of latitude (User-input)
 * \param[in] lat_end is end value of longtitude (User-input)
 * \param[in] lon_start is start value of longtitude (User-input)
 * \param[in] lon_end is end value of longtitude (User-input)
 * \param[in] time_start_index is start index of time (User-input)
 * \param[in] time_end_index is end index of time (User-input)
 * \param[in] time_step is time step (User-input)
 * \param[in] level is level will be processed (User-input)
 * \return void
 * \note this is a private helper function
 * \warning note that parameters are pointer
 */
static void set_domain_properties( float lat_start, float lat_end, float lon_start, float lon_end, int time_start_index, int time_end_index, int time_step,  int level) {
	//Check initial values
	checkLatLon(&lat_start, &lat_end, info.attr.lat.actual_range);
	checkLatLon(&lon_start, &lon_end, info.attr.lon.actual_range);
	checkTime(&time_start_index, &time_end_index,&time_step);
	if( level<0 || level>info.dims.level.len ) level = 0;
	//
	dom.lat.start = lat_start;
	dom.lat.end = lat_end;
	dom.lat.step = dom.dims.lats[0] - dom.dims.lats[1];/* Set the grid delta latitude */
	dom.lat.index.start = LAT_TO_IND(lat_end, dom.lat.step);/* Set the latitude start index for calculations */
	dom.lat.index.end = LAT_TO_IND(lat_start, dom.lat.step);/* Set the latitude end index for calculations */
	/*------------------------------------------------------------------------------------------------------------------------*/
	dom.lon.start = lon_start;
	dom.lon.end = lon_end;
	dom.lon.step = dom.dims.lons[1] - dom.dims.lons[0]; /* Set the grid delta longtitude */
	dom.lon.index.start = LON_TO_IND(lon_start, dom.lon.step); /* Set the longtitude start index for calculations */
	dom.lon.index.end = LON_TO_IND(lon_end, dom.lon.step);/* Set the longtitude end index for calculations */
	/*------------------------------------------------------------------------------------------------------------------------*/
	dom.time.index.start = time_start_index;
	dom.time.index.end = time_end_index;
	dom.time.step = time_step;
	/*------------------------------------------------------------------------------------------------------------------------*/
	dom.level = level;
}

/*!
 * \brief Prints the found cyclones to the standart output (file or screen)
 *
 * @param[in] cyc This is found cyclone list (double-linked)
 * @param[in] cyclone_type is the character array of cyclone type
 *
 * \return void
 * \note this is a private helper function
 * \warning note that parameters are pointer
 */
static void print_cyclones(struct CYCLONE* cyc, const char* cyclone_type) {

	printf(
	"\n%s Cyclones:\n"
	"\tRaw Count : %d\n"
	"\tEliminated : %d\n"
	"\tNet Count : %d {\n\n",
	cyclone_type,
	cyc->count.raw,
	cyc->count.eliminated,
	cyc->count.net);

	struct cyclone_item* p = cyc->first;
	register int i = 0;
	while(p) {
		i++;
		printf("%d\t%12.8g;%12.8g;%12.8g\n", i, p->lat, p->lon, p->height);
		p=p->next;
	};
	printf("}\n");
}

 /*!
 * \brief Calculates the accurate real min point by using the value of the grid point
 *
 * This function mainly gets a pointer to the grid point value (\a pGridPoint).
 * This value is a minimum candidate and had been compared against 8-neighbour point values around it before.
 * Using by \a scaleFactor and \a addOffeset in the function \ref hgt_find_real_min_location, real accurate minimum point is determined.
 *
 * \sa hgt_find_real_min_location
 * \sa cyclone_create
 *
 * @param[in] *pgridPoint This is the pointer to the value of the grid point.
 * @param[in] lp_lat This is the lattitude index of the grid point
 * @param[in] lp_lon This is the longtitude index of the grid point
 * @param[in] time This is the current time value. (Read from the NC file by time-index)
 *
 * \return cyclone_item* This is the new cyclone item created by the function
 * \note this is a private helper function
 * \warning note that *pgridPoint is a pointer
 */
static struct cyclone_item* get_accurate_min_point(short* pgridPoint, int lp_lat, int lp_lon, double time ) {
	double x_loc = 0, y_loc = 0, val = 0;
	struct cyclone_item* new_cyc = NULL;
	log_deb("\nlat = %.1f\tlon=%.1f{",IND_TO_LAT(lp_lat),IND_TO_LON(lp_lon));
	/* Now we need to calculate the real min point and value at the moment at this point */
	if( hgt_find_real_min_location(pgridPoint, info.attr.hgt.scaleFactor, info.attr.hgt.addOffset, &x_loc, &y_loc, &val) == 0 ) {

		/* Caculated x_loc and y_loc values are relative to mid-point, so we had to convert them to global coordinate values */
		double lat_new = IND_TO_LAT(lp_lat) + y_loc;
		double lon_new = IND_TO_LON(lp_lon) + x_loc;

		new_cyc = cyclone_create(time, dom.level, val, lat_new, lon_new);
	}
	log_deb("}");
	return new_cyc;
}

/*!
 * \brief Compares the new cyclone to other cyclones found before and eliminates the lower one
 *
 * This Function compares the new cyclone against the other cyclones in the \ref cyc "double-linked list(cyc)" found before and
 * removes shallow ones around
 *
 * \sa Haversine_Distance
 * \sa cyclone_overwrite
 * \sa cyclone_delete
 *
 * @param *new_cyc is the new cyclone item just found
 * @param[out] *cycs is the cycones that found before.
 * @param distance_limit is the desicion paramater to eliminate the cyclone or not.
 *
 * \return an integer value that a close cycone to the new one had been found or not.
 * \retval '0' - No close cyclone found
 * \retval '1' - Another cyclone found
 * \note this is a private helper function
 * \warning note that \a *new_cyc and \a *cycs is a pointer
 */
static int eliminate_close_cyclones(struct cyclone_item* new_cyc, struct CYCLONE *cycs, const double distance_limit) {
	struct cyclone_item* p = cycs->first;
	struct cyclone_item** first = &(cycs->first);
	struct cyclone_item** last = &(cycs->last);
	int isClose = 0; // Means no close cyclone found
	while(p) { /* Compare the new cyclone to others */
		/* Calculate the distance between new cyclone and others on the same level */
		double dis = Haversine_Distance(p->lat, p->lon, new_cyc->lat, new_cyc->lon );
		if(dis <= distance_limit) {/* If distance is smaller then or equal distance_limit (km.) */
			/* If new_cyclone height is lower than one of cyclones */
			if( (new_cyc->height) < (p->height) ) {
				(cycs->count.eliminated)++;/* Increase eliminated cyclone count */
				if(isClose==0) {/* First close cyclone determined */
					cyclone_overwrite(new_cyc, &p, first, last);/* Swap the new cyclone with other higher one */
					isClose = 1;/* We found first cyclone point close to new one. */
				}else {/* Another close cyclone determined - REMOVE IT! */
					cyclone_delete(p, first, last);
				}
			}
		}
		p = p->next;
	}
	if(isClose==0) {
		/* If there is no cyclone close together, then add the new found cyclone to the List */
		if(*first==NULL) *first = new_cyc;
		cyclone_add(new_cyc, last);
	}
	return isClose;
}

/*!
 * @brief Process the 3D data array for a given specific time and searchs for the cyclone points.
 *
 * This Function scans a 3D data array for a given spcific \a time. To scan, checks if every
 * grid points are minimum candidate or not by \ref hgt_is_min function. If a minimum candidate is found,
 * then calculates the accurate minimum coordinate by \ref get_accurate_min_point function.
 * After that, new cyclone is compared against other cyclones on the same chart for elimination by
 * \ref eliminate_close_cyclones function and prints them to the selected output (screen or file)
 *
 * @sa hgt_is_min
 * @sa get_accurate_min_point
 * @sa eliminate_close_cyclones
 * @sa print_cyclones
 *
 * @param ***data has 3D geoptantial height data for given time from the file.
 * @param time is the time step for read process
 *
 * @return void
 * @note this is a private helper function
 * @warning note that \a ***data is pointer of pointer of pointer. (Represents 3D array)
 */
static void scan_level(short*** data, double time) {

	/* Variables *///{

	register int lp_lat, lp_lon;
	int lat_start = dom.lat.index.start;
	int lat_end = dom.lat.index.end;
	int lon_start = dom.lon.index.start;
	int lon_end = dom.lon.index.end;
	int level = dom.level;

	struct CYCLONE chart_cyc;
	chart_cyc.first = chart_cyc.last = NULL;
	chart_cyc.count.raw = 0;
	chart_cyc.count.eliminated = 0;
	chart_cyc.count.net = 0;
	struct COUNT* count = &(chart_cyc.count);
	//}

	/*------------------------------------------------------------------------------------------------------------------------*/
	loop(lp_lat, lat_start, lat_end) {
		loop(lp_lon, lon_start, lon_end) {
			short* pgridPoint = &data[level][lp_lat][lp_lon];
			if( hgt_is_min(pgridPoint) ) {
				/* We need to calculate the real min point and value at the moment at this point*/
				struct cyclone_item* new_cyc = get_accurate_min_point(pgridPoint, lp_lat, lp_lon, time);
				if(new_cyc != NULL) {
					count->raw++;
					eliminate_close_cyclones(new_cyc, &chart_cyc, ELIMINATION_RADIUS);
				}
			}
		}
	}
	/*------------------------------------------------------------------------------------------------------------------------*/

	count->net = count->raw - count->eliminated;
	connect_dls(&(cyc.first), &(cyc.last), &(chart_cyc.first), &(chart_cyc.last));
	cyc.count.raw += count->raw;
	cyc.count.eliminated += count->eliminated;
	cyc.count.net += count->net;

	print_cyclones(&chart_cyc,"Chart");
}

/*!
 * @brief Sets initial conditions and starts the search process for the cyclones
 *
 * This Function sets initial conditions given by the user input and creates 3D data array
 * to read the geopotantial height values by the \ref createDataArray function. Then
 * loops in given time range by time-step and reads domain data from file for specified
 * level by \ref NC_get_data function. Scans the grid points in the 3D data array and
 * searches for minimum grid points by \ref scan_level function.
 *
 * @sa NC_get_data
 * @sa scan_level
 * @sa print_cyclones
 *
 * @return void
 * @note this is a private helper function
 */
static void calculate() {

	/* Set Initial Conditions -----------*///{
	int time_start = dom.time.index.start;
	int time_end = dom.time.index.end;
	/* The start and count arrays will tell the netCDF library where to read our data. */
	size_t* start = arr_malloc(size_t, info.dims.count);
	size_t* count = arr_malloc(size_t, info.dims.count);
	/* Read the data. Since we know the contents of the file we know
    * that the data arrays in this program are the correct size to
    * hold one timestep. */
	count[0] = 1;
	count[1] = info.dims.level.len;
	count[2] = info.dims.lat.len;
	count[3] = info.dims.lon.len;
	start[1] = 0;
	start[2] = 0;
	start[3] = 0;
	/* Program variables to hold the data we will read. We will only
	need enough space to hold one timestep of data; one record. */
	short*** hgt_in = createDataArray(info.dims.level.len, info.dims.lat.len, info.dims.lon.len);
	//}

	register int lp_time;/* Loop index for time. */
	loop(lp_time,time_start,time_end) {
		start[0] = lp_time;
		double time = dom.dims.times[lp_time];
		printf("\nTime = %d {\n",lp_time);
		log_deb("Time: %d {",lp_time);
		NC_get_data(hgt_in[0][0], start, count);
		scan_level(hgt_in, time);
		log_deb("\n}\n");
		printf("\n}\n");
	}
	//
	print_cyclones(&cyc, "Total");

	free(start);
	free(count);
	free(hgt_in);
}

//}

/* Public Routines ----------------------------------------------------------------------------------*///{

/** @name Public
 *  These functions are defined in the header file \ref Siklon.h
 */
//@{

/*!
 * @brief Initialize the needed variables and resources to start the siklon search process.
 *
 * Opens the specified NC file, read and checks the file properties
 * against user input and sets the domain properties. This function must
 * be used first of all.
 *
 * @sa NC_open_file
 * @sa NC_get_variable_ids
 * @sa NC_read_dimension_lengths
 * @sa NC_read_attributes
 * @sa NC_read_coord_var_data
 * @sa set_domain_properties
 *
 * @param *fileName is NC file to be read.
 * @param lat_start is given start index for the latitude
 * @param lat_end is given end index for the latitud
 * @param lon_start is given start index for longtitude
 * @param lon_end is given end index for longtitude
 * @param time_start is given start vindex value for time range
 * @param time_end is given end index value for the time range
 * @param time_step is the time step
 * @param level is the level will be processed
 *
 * @return if no error, returns 1.
 */
int sk_Init(char *fileName, float lat_start, float lat_end, float lon_start, float lon_end,
			int time_start, int time_end,int time_step, int level) {
	NC_open_file(fileName);
	NC_get_variable_ids();
	NC_read_dimension_lengths();
	NC_read_attributes();
	NC_read_coord_var_data();
	set_domain_properties(lat_start, lat_end, lon_start, lon_end, time_start, time_end, time_step, level);
	cyc.first = cyc.last = NULL;
	return 1;
}

/*!
 * @brief Starts the siklon search process.
 *
 * @sa calculate
 *
 * @return if no error, returns 1.
 * @warning This public function must be used after the \ref sk_Init function.
 */
int sk_scan(){
	calculate();
	return 1;
}

/*!
 * @brief Releases all allocated resources.
 *
 * Finishes the search process and releases all used resources.
 *
 * @sa calculate
 *
 * @return if no error, returns 1.
 * @warning This public function must be used after the \ref sk_scan function.
 */
void sk_finalize(){
	nc_close(info.id.nc);
	struct cyclone_item* p = cyc.first;
	struct cyclone_item* next = NULL;
	while(p) {
		next = p->next;
		free(next);
		p = next;
	}
}

/*!
 * @brief Prints the NC file info to the standart output
 *
 * @param *fileName is the NC file will be processed
 *
 * @sa NC_open_file
 * @sa NC_get_variable_ids
 * @sa NC_read_dimension_lengths
 * @sa NC_read_attributes
 * @sa NC_read_coord_var_data
 * @sa NC_print_coord_var_data
 *
 * @return if no error, returns 1.
 * @warning This public function must be used after the \ref sk_scan function.
 */
void sk_print_file_info(char* fileName) {
	NC_open_file(fileName);
	NC_get_variable_ids();
	NC_read_dimension_lengths();
	NC_read_attributes();

	/** Print File Header Info *///{
	(void)fprintf(stdout,
	"open_NC_file:\n"
	"\tid.nc = %d \n"
	"get_variable_ids:\n"
	"\tid.level = %d\n"
	"\tid.lat = %d\n"
	"\tid.lon = %d\n"
	"\tid.time = %d\n"
	"\tid.hgt = %d\n"
	"file_Header_Info\n"
	"\tdims: Count = %d\n"
	"\t  lon:\tid=%d  len=%zd\n"
	"\t  lat:\tid=%d  len=%zd\n"
	"\t  level:id=%d  len=%zd\n"
	"\t  time:\tid=%d  len=%zd\n"
	"\tattributes:\n"
	"\t  level:\n"
	"\t    units : \"%s\"\n"
	"\t    actual_range: %f, %f\n"
	"\t    long_name: \"%s\"\n"
	"\t    positive: \"%s\"\n"
	"\t    GRIB_id : %d\n"
	"\t    GRIB_name : \"%s\"\n"
	"\t  lat:\n"
	"\t    units : \"%s\"\n"
	"\t    actual_range: %f, %f\n"
	"\t    long_name: \"%s\"\n"
	"\t  lon:\n"
	"\t    units : \"%s\"\n"
	"\t    actual_range: %f, %f\n"
	"\t    long_name: \"%s\"\n"
	"\t  time:\n"
	"\t    long_name: \"%s\"\n"
	"\t    units : \"%s\"\n"
	"\t    delta_t: \"%s\"\n"
	"\t    actual_range: %f, %f\n"
	"\t  hgt:\n"
	"\t    long_name: \"%s\"\n"
	"\t    actual_range: %f, %f\n"
	"\t    valid_range: %f, %f\n"
	"\t    units : \"%s\"\n"
	"\t    addOffset : %f\n"
	"\t    scaleFactor : %f\n"
	"\t    missing_value : %d\n"
	"\t    GRIB_id : %d\n"
	"\t    GRIB_name : \"%s\"\n"
	"\t    var_desc : \"%s\"\n",
	info.id.nc,
	info.id.level, info.id.lat, info.id.lon, info.id.time, info.id.hgt,
	info.dims.count,
	info.dims.lon.id, info.dims.lon.len,
	info.dims.lat.id, info.dims.lat.len,
	info.dims.level.id, info.dims.level.len,
	info.dims.time.id, info.dims.time.len,
	info.attr.level.units,
	info.attr.level.actual_range.min,info.attr.level.actual_range.max,
	info.attr.level.long_name,
	info.attr.level.positive,
	info.attr.level.GRIB_id,
	info.attr.level.GRIB_name,
	info.attr.lat.units,
	info.attr.lat.actual_range.min, info.attr.lat.actual_range.max,
	info.attr.lat.long_name,
	info.attr.lon.units,
	info.attr.lon.actual_range.min, info.attr.lon.actual_range.max,
	info.attr.lon.long_name,
	info.attr.time.long_name,
	info.attr.time.units,
	info.attr.time.delta_t,
	info.attr.time.actual_range.min, info.attr.time.actual_range.max,
	info.attr.hgt.long_name,
	info.attr.hgt.actual_range.min, info.attr.hgt.actual_range.max,
	info.attr.hgt.valid_range.min, info.attr.hgt.valid_range.max,
	info.attr.hgt.units,
	info.attr.hgt.addOffset,
	info.attr.hgt.scaleFactor,
	info.attr.hgt.missing_value,
	info.attr.hgt.GRIB_id,
	info.attr.hgt.GRIB_name,
	info.attr.hgt.var_desc);
	//}

	NC_read_coord_var_data();
	NC_print_coord_var_data();

	nc_close(info.id.nc);
}

/*!
 * @brief Prints the initial conditions set by user in the command-line paramters
 * @return void
 */
void sk_printInitialConditions(){
	(void)fprintf(stdout,
	"Initial Conditions\n"
	"lat:\n"
	"  Start:%f\tIndex:%d\n"
	"  End :%f\tIndex:%d\n"
	"lon:\n"
	"  Start:%f\tIndex:%d\n"
	"  End :%f\tIndex:%d\n"
	"time:\n"
	"  StartIndex:%d\n"
	"  EndIndex :%d\n"
	"  Step :%d\n"
	"level: %d\n\n",
	dom.lat.start, dom.lat.index.start, dom.lat.end, dom.lat.index.end,
	dom.lon.start, dom.lon.index.start, dom.lon.end, dom.lon.index.end,
	dom.time.index.start, dom.time.index.end, dom.time.step,
	dom.level);
}

//@}

//}











